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Abstract. In the linear collisionless limit, a zonal potential perturbation in a toroidal 
plasma relaxes, in general, to a non-zero residual value. Expressions for the residual 
value in tokamak and stellarator geometries, and for arbitrary wavelengths, are derived. 
These expressions involve averages over the lowest order particle trajectories, that 
typically cannot be evaluated analytically. In this work, an efficient numerical method 
for the evaluation of such expressions is reported. It is shown that this method is 
faster than direct gyrokinetic simulations performed with the Gene and EUTERPE 
codes. Calculations of the residual value in stellarators are provided for much shorter 
wavelengths than previously available in the literature. Electrons must be treated 
kinetically in stellarators because, unlike in tokamaks, kinetic electrons modify the 
residual value even at long wavelengths. This effect, that had already been predicted 
theoretically, is confirmed by gyrokinetic simulations. 
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1. Introduction 

Whereas the reduction of turbulent transport by zonal flows is a widely accepted 
phenomenon [1, 2], its quantitative understanding is still poor. The determination 
of the zonal flow amplitude and turbulent transport level are nonlinear questions whose 
answers require costly gyrokinetic simulations. The cost grows enormously if, by means 
of parameter scans, one wants to know how those quantities depend on the magnetic 
geometry or the plasma conditions. A useful simplihcation is provided by the initial 
value problem consisting of calculating the long-time and collisionless evolution of a 
zonal perturbation. Among other reasons, its usefulness is due to the fact that an exact 
expression for the value of the perturbation at f = oo, called residual value, can be 
derived. Even though the explicit evaluation of the final expression can only be carried 
out in simplihed geometries and for long wavelengths of the perturbation, the analytical 
solution gives insight into the physics. This partly explains the attention attracted by 
the work of Rosenbluth and Hinton [3] and the effort put on subsequent extensions 
that we cite below. In stellarator research, the interest in this problem was spurred by 
the suggestion [4] of a direct relation between zonal flow residual value and turbulent 
transport level. 

The seminal paper [3] dealt with long-wavelength potential perturbations in 
large aspect ratio and circular cross section tokamaks. Explicit solutions of related 
problems for more complex tokamak geometries and arbitrary wavelengths have been 
given in [5, 6]. In recent years, several articles have addressed the problem in 
stellarators [7, 8, 9, 10, 11]. In this paper we report on a code to evaluate fast 
and accurately the exact expressions of the zonal flow residual value, for arbitrary 
wavelengths and for tokamak and stellarator geometries. But why is this useful if the 
same answer can be obtained, in principle, by linear runs of a local gyrokinetic code? 
As we will illustrate later on, it turns out that the solution by means of gyrokinetic 
simulations of the residual zonal flow problem, especially at short wavelengths, is very 
demanding in terms of computational resources. We will show that our method is faster. 
These differences can be of several orders of magnitude in computing time, especially 
in the case of stellarators. Moreover, the code and the results of this work are not only 
interesting from the point of view of physics, but also for validation of gyrokinetic codes. 

The rest of the paper is organized as follows. In Section 2 we derive the residual 
zonal flow expressions for arbitrary wavelengths and magnetic geometry. In Section 3 
we report on the code used to evaluate the expressions presented in Section 2. As 
a check, we compare our calculations with analytical results from [5, 6], obtained in 
simplihed tokamak geometry. In order to avoid any confusion, we note that the results 
in [5, 6] were not derived as an initial value problem, but as the stationary solution of 
a forced system. We explain this in more detail in Section 3. In Section 4 our results 
are compared with local and global gyrokinetic simulations of the initial value problem, 
employing the codes Gene [12, 13, 14, 15] and EUTERPE [16, 17], respectively. We 
also compare the differences in computational time required by each approach showing 
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that the method introduced in this paper is faster than those gyrokinetic simulations. 
Stellarator residual values are calculated, using our method and also gyrokinetic codes, 
for a range of wavelengths much wider than previously available in the literature. The 
stellarator calculations are done for the standard conhguration of Wendelstein 7-X (W7- 
X). We comment on a purely stellarator effect already predicted in [ 8 , 10]; namely, that 
the approximation of adiabatic electrons is always incorrect (even at long wavelengths) 
for the purpose of determining the residual zonal flow in stellarators. We quantify the 
error by computing, with our method, the residual value when kinetic or adiabatic 
electrons are used. This result is conhrmed by gyrokinetic simulations. The conclusions 
are presented in Section 5. 

Throughout this paper, we assume that the electrostatic potential ip associated to 
the zonal flow is constant on flux surfaces. It is true that zonal flows {i.e. flows that 
are only weakly damped, and therefore remain in the plasma for long times) usually 
correspond to electrostatic potentials with small variations on flux surfaces, but we 
should emphasize that assuming that ip is constant on flux surfaces is, at most, a good 
approximation. 

2. Linear collisionless evolution of zonal flows 

In this section, we give a detailed calculation of the residual zonal flow for arbitrary 
wavelengths in tokamak and stellarator geometries. We solve the linear and 
collisionless gyrokinetic equations at long times, assuming that the electrostatic potential 
perturbation is constant on flux surfaces. 

In strongly magnetized plasmas, one employs the smallness of = pts/L to 
average over the gyromotion. Here, L is the characteristic length of variation of the 
magnitude of the magnetic held B, pts = Vts/klg is the thermal gyroradius, Vts = \jTs/ms 
is the thermal speed, = ZgeB/mg is the gyrofrequency, Tg is the equilibrium 
temperature, rUg is the mass, and ZgC is the charge of species s, where e is the proton 
charge. Gyrokinetic theory [18, 19, 20, 21 , 22 , 23] gives a procedure to rigorously derive 
the gyroaveraged kinetic equations order by order in <C 1. The averaging operation 
is conveniently expressed in a new set of phase space coordinates, called gyrokinetic 
coordinates. Denote by {r, v} the particle position and velocity. The coordinate 
transformation, to lowest order in pt**, is given by 

r = R + p,(R, n. A, 7 ) + 0{pl^L), 

V = v\\ (R, V, A, cr)b(R) + klgPg{R, v, A, 7 ) x b(R) + 0{ptgi,vts). (1) 

In equation (1), R is the gyrocenter position, v is the magnitude of v, A = B~^v\/v‘^ is 
the pitch angle, a = n||/|n||| is the sign of the parallel velocity 

ny (R, n. A, a) = av ^/l-XB(R), (2) 

b is the unit vector in the direction of the magnetic held B, is the component of the 
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TTl V \ 

= ^J-^^[e 2 (R)cos 7 -ei(R)sin 7 ]. (3) 

Here, ei(R) and 62 (R) are unit vector fields orthogonal to each other which satisfy 
ei X e 2 = b at every point. Finally, the gyrophase 7 is 

7 = arctan(v ■ e 2 /v ■ ei). (4) 

We introduce straight field line coordinates {'0,6*,^}, where ^|J G [0,1] is the radial 
coordinate defined as the normalized toroidal flux -0 = 6 ^ is a poloidal angle 

and C is a toroidal angle, with 0,(^6 [0,1). The magnetic field in these coordinates is 
written as 


B = -vl/;(0)V0x V(C-g(0)0). (5) 

Here, qigp) = 4/^('0)/4/p(0) is the safety factor, and and 4/^(0) are the derivatives 

of the poloidal and toroidal fluxes with respect to 0. It will be convenient to define the 
coordinate a \= C, — that labels magnetic field lines on each flux surface. Unless 

otherwise stated, we use {0, 9, a} as the set of independent spatial coordinates. Note 
that B-V0 = O, B-Va = 0. We employ 6 as the coordinate along a field line. 

The distribution function in gyrokinetic coordinates, Fg = Fs{'ip,9,a,v, Xyay'jfl), 
can be written as 


F0R,n, A,cr, 7 ,t) = Fso(R,v) + F,i(R, n. A, cr, f) + 0{p^g^Fso), (6) 

where Fgi = 0{pts*Fso) and Fgo is a Maxwellian distribution whose density Ug and 
temperature Tg = mgV^g are flux functions. 


Fgo(R,v) := 


n00(R)) 


exp 


(7) 


(x/^nt00(R)))3 V 2n200(R)) 

The lowest-order quasineutrality condition implies = 0. Note that to 

0 (pts*Fso) the distribution function is independent of the gyrophase. 

The linear and collisionless time evolution of Fgi is given by [24, 23] 

dtHgi + (nil b + Vd0 ■VHgi = -^dt{(p)Fgo 


+— ( V(V7) X b) ■ V0 


— + 
Ug 


nigV 


2 T. 


sOj 


( 8 ) 


where primes denote differentiation with respect to 0 , the function Hgi is defined by 


Hsi — Fgi + -—{p)Fgo, 

S 

the gyroaveraged electrostatic potential is 

1 

(ip) (R, V, ~ ^ / V' (R + ft(R. f. 7 ), t) d7 


(9) 


( 10 ) 
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and the magnetic drift velocity reads 


vwo = —b X 


(1 - \B)h ■ Vh +-VB 


( 11 ) 


The orderings in gyrokinetic theory allow to separate the variations of the fields 
on the small and large scales, and decompose in Fourier modes with respect to the 
former. Since we are interested in studying the evolution of an electrostatic potential 
perturbation that depends only on and the problem is linear, we can take a single 
mode of the form 


(p{r,t) = ipk{'il^{r),t)exp{ik^'il:{r)). 

Here, 


( 12 ) 

(13) 


with /d_l(R) = fcpiV'^(R)|, and varies on the macroscopic scale L. Observe that, due 
to the effects of magnetic geometry, the dependence of on cannot be avoided even 
for flat density and temperature profiles. A recent explanation of scale separation, as 
well as a proof of the equivalence between the local and global approaches to gyrokinetic 
theory can be found in [25]. 

To lowest order, the gyroaveraged electrostatic potential is 


((p)(R,n, A,f) = (pfc(V’(R),t) Jo{k±ps) exp{ik^i/j{R)), 
where the magnitude of the gyroradius vector is 

rUsV I A 


PsiR^v, A) = 


ZgC y R(R) 


and Jo is the zeroth-order Bessel function of the first kind. 


(14) 


(15) 


1 

Jo{x) = — / exp(ia; sin 7 )d 7 . (16) 

27r Jo 

If the electrostatic potential has the form (12), then the distribution function can be 
written as 


FsiiR, V, A, a, t) = fs(pp(R), 9(R), q((R), n. A, a, t) exp(i/c^yi(R)) (17) 

and consequently 


Jf^i(R,n, A,cr, t) 


hs{'ip(R),9(R), a(R), n. A, a, t) exp(i/c^V>(R)), 


(18) 


where fs and hg vary on the scale L. Then, equation (8) becomes 
T n|| b ■ V T ik-ipui^ hg dtPkJosFso^ 


(19) 


where we have used the notation Ug := ■ 'Vfl> for the radial magnetic drift frequency 

and Jqs = Jo{k±pg). From now on, and for brevity, we omit the dependence of ipk on -0; 
that is, we write (pk(t) instead of ipk{fl’,t). 

Denote by u the frequency associated to the time derivative in (19). The objective 
is to expand (19) in powers of oj/{ytsL~^) <C 1, solve the lowest order equations and 
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determine (pflt) in the limit t —)■ oo. The u/{vtsL~^) ^ 1 expansion means, in particular, 
that we average over the lowest order particle trajectories and solve for time scales much 
longer than a typical orbit time, which is 0{L/vts)- We dehne the orbit average for a 
phase-space function Q(V’, 0 ,«, v, A, a, t) as 

{B Q/\v\\\)^/{B/\v\\\)^ for passing particles 

ojb f ddQ/(vii b ■ V^) for trapped particles, 

where uj/j := [^d6'/(u|| b ■ is the bounce frequency. Given a function Gib'll!, 9, a), 

the flux surface average is dehned by 


Q:= 


( 20 ) 


(G), = vy,) 


-1 


( 21 ) 


d6 day/g G{i/j,9,a). 

Jo Jo 

Here, y/g = [(VV' x V^) ■ is the square root of the metric determinant and 

= Jq d^ Jq da-s/g is the derivative of the volume enclosed by the flux surface 
labeled by The symbol f stands for integration over the trapped trajectory, where 
the bounce points df, are the solutions of 1 — \B{'ip, Of,, a) = 0 for given values of and 
a, and given an initial condition for the particle trajectory. 

Observe that the orbit average operation has the property 


U||b-VQ = 0 (22) 

for any single-valued function Q. We write the radial magnetic drift frequency as a sum 
of its orbit averaged and fluctuating parts 

Co’s = -|- U|| b ■ V5s, (23) 

where 6s = Sgigp, 9, a, v, A, a), that we choose to be odd in a, is the radial displacement 
of the particle’s gyrocenter from its mean flux surface. The solution of the magnetic 
differential equation that determines 6s is given in Appendix A. 

Dehning := hs exp{ik^6s) and ipk '■= (fk exp {iky 6 s), equation (19) yields 


-F U||b ■ V -f ikyUJs'^ hs = ^dt(£kJosFso- 
It is worth noting that the expansion in u/{vtsL~^) only makes sense if 

kyUJ^ OJ 


For kypts 


VtsL-^ VtsL-^ 

1 , this implies 

6J~s OJ 

OJ. 


< 1 . 


< 1 . 


(24) 


(25) 


(26) 


VuL~' 

This trivially holds in a tokamak because = 0 for all trajectories. In a generic 
stellarator, = 0 only for passing particles. Then, condition (26) requires that the 
secular radial drifts of trapped particles be sufficiently small. We assume that this is 
the case and carry out the expansion in oj/{vtsL~^). 

We write 


h, = Af + AW + Af > + ..., 


(27) 
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with ~ (jjflvtsL ^). Then, we expand equation (24). To lowest order, one 

obtains 




b ■ Vh® = 0, 


implying that is constant along the lowest order trajectories; i.e. 


lgl=h^. 

To next order, we have 




{dt + ik^Us) + nyb ■ = -^dtcpkJosF, 

-L .« 


sO- 


(28) 


(29) 


(30) 


We do not write 93 ® to ease the notation. The orbit average of (30) annihilates the 
term nyb ■ V/^^, and we find 


{dt + ik^iVs) = ^dt^JosFso. 


(31) 


It is useful to work in Laplace space in order to solve this equation. The Laplace 
transform of a function Q{t) is defined as Q{p) = Q{t)e~P^dt, where p denotes the 
variable in Laplace space. We apply it to (31) and obtain 


Z 6 

(p + ik^uJ^) = -j^pflkJosFso + /s(0). 

J- ft 


(32) 


Here, /s(0) := fs{0) exp{ik^6s) and /s(0) is the initial condition for fg] i.e. /s(0) = 
/s('0,6',Q;,n, A,a,0). 

The solution of (32) yields 


0 






(33) 


In order to have a closed system of equations we employ the gyrokinetic 
quasineutrality equation (see, for example, references [24, 23]), 




= 




(34) 


Here, the short-hand notation J Qd^v means, for a function Q('0, 9, a, v, A, a, 7 ), 


/»27r 

Qd\.= W / d7 

<x=-W0 


v'^B 


Using (12) and (18), and flux-surface averaging, we get 


E 


Tg 


‘^s Pk 


^z.jj„.h.d= 


(35) 


(36) 
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To lowest order in u/{vtsL -C 1, and after transforming to Laplace space, equation 
(36) gives 


E 




'^s — 




(37) 


We employ (33) to write the right side of (37) in terms of the electrostatic potential and 
the initial condition, and solve for ipk. The result is 


^fc(p) = 






p+ik^uj. 




Os 


(38) 


where we have simplified the notation by defining 


/ 1 roo p 

{Q}s - (E / / 

\^=-lJo Jo 


dA -^===Q{'il!, 9, a, V, A, a) F^o 


(39) 




for gyrophase independent functions on phase space. 

The residual value is found from the well-known property of the Laplace transform 


lim (pk{t) = limpflflp)- 

t^oo p^O 

Applying (40) to equation (38), we find 


E.z.- 

0 ^k^Ss 

JosF^^^^fs{0)/F,o] 

LOs=0 

S 

Z^s Ts 



\ lJ^=0 

J s 

e Jqs 


(40) 


(41) 


where (pk{oo) = \imt^ooPk(t). The superscript = 0 means that the integration 
is performed only for particles whose trajectory satisfies aJj" = 0. In tokamaks, this 
property holds true for both trapped and passing particles, and therefore the integrals 
in (41) are performed over the whole phase space. In stellarators, = 0 is satisfied 
exclusively for passing particles. Only in perfectly omnigenous stellarators [26, 27, 28] 
have trapped particles vanishing average radial magnetic drift. Hence, in a generic 
stellarator, the integrals in (41) with superscript = 0 are performed only over the 
passing region of phase space. 

The residual level is usually dehned as the normalized value ipk{oo)/ipk{0). The 
relation between /s(0) and (pfc(O) is given by the flux-surface averaged quasineutrality 
equation at t = 0 , 

^^n,(l-r„(fcipf,))^^»(0) = ^^z,|j„./,(0)d=>t.^ , (42) 

Here, we have employed the identity J Jg (fc_LPs)Tosd^n = nsTo^k^pffl, where 
^oiklpi) := loiklpl), and Jg is the zeroth order modified Bessel function. We 

will use the notation Tg^ = To^k'^pffl. 
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For later comparison with gyrokinetic simulations, it will be useful to have at hand 
the expressions corresponding to the approximation of adiabatic electrons. Using this 
approximation, equation (36) can be written as 


E 


—Us^kip) = 



(43) 


Proceeding as shown previously for the fully kinetic case, we find that the expression 
for (pk{oo) reads 


Pk{oo) 


z.. 

^—ik■^|J6s 

Jo,e'"v-5^/,(0)/F,o} 

UJs=0 

s 

Ts 



\ ij^=0 

) s 

e Jqs 


(44) 


In this case, the equation that relates /s(0) and (pk{0) is 

'^^ns{l-Tos)^Pk{0) = f Josfs{0)dh) ■ (45) 

s^e ^ \ SJ^e / ^ 

The residual zonal flow was computed in reference [3] for k^pu 1 in a large 
aspect ratio, circular cross section tokamak with adiabatic electrons. An extension of 
the derivation of [3] was proposed in references [5, 6] to allow for short-wavelength 
perturbations, and also for kinetic electrons and more complex tokamak geometries. 
In reference [6], comparisons of analytical calculations with gyrokinetic simulations are 
shown. The enhancement of tokamak residual zonal flows at short wavelengths was 
originally found in reference [12] by means of gyrokinetic simulations with the code 
GS2. But the short-wavelength calculations of [12, 5, 6] do not correspond to an initial 
zonal value problem because the quasineutrality equation is forced with a source term. 
At long wavelengths, the initial value problem and the forced system give the same 
result. This will be explained in more detail in Section 3. 

In stellarators, the residual zonal flow calculation has been done in [7, 8, 9, 10, 11]. 
In references [9, 10, 11] the emphasis is put on long-wavelength zonal flows. In references 
[7, 8], the derivation of the equations is valid for long and short-wavelengths but some 
approximations are used to describe the magnetic geometry. 


3. Evaluation of the expressions for the residual zonal flow 

The evaluation of the right side of equation (41) (and, of course, (44)) requires the 
calculation of quantities of the form 

' Ip 

for functions P = P(7, v, A, a) and Q = Qigp, 9, a, v, A, a), where the orbit average 
is defined in (20). These averages depend on the details of the magnetic field and cannot 
be evaluated analytically, except in simplified cases (for example, in the cases considered 
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in [3, 5, 6]). In this work, we evaluate equation (41) using the framework of the code 
CAS3D-K [29, 30]. For this purpose, we have included in this code the relevant hnite 
Larmor radius effects, the solution to the magnetic differential equations described in 
Appendix A, and the integration over the velocity coordinates v and a. 

The code CAS3D-K is well suited to perform the average (20). Since the lowest 
order particle trajectories lie entirely on flux surfaces, all the calculations are local, 
thus permitting a parallelization by flux surface using MPI. The magnetic equilibrium 
is obtained from the 3D MHD equilibrium code VMEC [31] and then transformed 
into Boozer coordinates which can easily be transformed to the coordinates 

On a given flux surface, the pitch angle A distinguishes between passing 
and trapped trajectories. The passing-trapped boundary is given by Ac = 
where is the maximum of B on the flux surface. Passing particles have A values 

with 0 < A < Ac and trapped particles are those with Ac < A < l/i?™“, where 
^mm jg minimum of B on the flux surface. Trapped particles can live inside 
one or several magnetic held periods. In CAS3D-K, they are grouped by the number 
of periods they go through. The groups are obtained by setting a large number of 
initial conditions for the trajectories, and hnding the bounce points from the bounce 
condition 1 —Ai?('0, 6*6, a) = 0 for constant xjj and a. From this procedure, the boundaries 
of each group are found and the numerical integration for a given group is performed 
by covering the region they dehne with new trajectories. Note that each group requires 
different numerical resolution. Some trapped trajectories close to the passing-trapped 
boundary may require a large number of periods until the bounce points are found. 
If this number is sufficiently large, typically larger than 500 periods, the trajectories 
are then considered as passing. A Gauss-Legendre quadrature scheme is used for the 
integration in 6*, a and A, which avoids the numerical problems that may appear at 
points where 1 — \B = 0. 

In reference [9], the phase space integrations with CAS3D-K discussed in the 
previous paragraph were employed to obtain the zonal flow frequency in stellarator 
geometry in the long-wavelength limit. In this limit, there are some simplihcations. We 
have implemented new functionalities in CAS3D-K that also allow to deal with short 
wavelengths. Next, we describe the main features of this extension of CAS3D-K. 

First, we note that the integration of the resulting expressions over the velocity 
coordinate v could be performed analytically in [9]. If one wants to calculate the averages 
involved in (41), the integration over the velocity coordinates v and a must be computed 
numerically. We have included in CAS3D-K the integration over these coordinates on 
top of the 6*, a and A integrations. For the integration over n, a linear scheme has 
been used. This scheme allows the computation of any moment in v of the Maxwellian 
distribution function. 

Second, the equations to obtain the residual level (41) and (44) incorporate hnite 
Larmor radius effects. These effects are encoded in Joik^ps), and exp(±iA;^(5s). 

Here, = kp\V^l^\ where kp is an input parameter and the quantity iV'i/'l = |V'0|(6*,Q!) 
is obtained from the VMEC equilibrium. We have included those factors and adapted the 
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Figure 1. Radial dependence of the residual level given by (47), by (48), and by the 
evaluation of (44) with CAS3D-K in the long-wavelength limit. A tokamak with major 
radius R = 1.7 m, minor radius a = 0.4 m, and q profile given in figure 2 has been 
used. 

resolution of the integration on phase space to their strong oscillatory behavior at short 
wavelengths. We have also implemented in the code the expression of 6s in stellarator 
geometry, for both passing and trapped particles (see the derivations in Appendix A). 
For passing particles, Ss is given by equation (A. 7) and for trapped particles it is given 
by equation (A. 17). 

The above modihcations included in CAS3D-K allow us to calculate the residual level 
in tokamak or stellarator geometry for arbitrary wavelengths. In the rest of the paper, 
we use the terminology “CAS3D-K” and “extension to CAS3D-K” interchangeably. 

As a preliminary check of these extensions to CAS3D-K, in this section we compare 
our results with analytical results available in the literature. For these comparisons, we 
take a plasma consisting of singly charged ions and electrons, and assume flat density 
and temperature prohles with the same values for both species. 

In reference [3], Rosenbluth and Hinton (R-H) calculated the residual level in large 
aspect ratio tokamaks with circular cross section and adiabatic electrons, in the limit 
k±pti ^ 1. Denote the safety factor by q and the inverse aspect ratio by e = {a/R)y/^, 
where a is the minor radius and R is the major radius. The result obtained in [3] is 

P>k{oo) _ 1 

(^fc(O) 1-b 1.6q'2£-V2 • 



( 47 ) 










Residual zonal flows in tokamaks and stellarators at arbitrary wavelengths 


12 


2.6 

2.4 
2.2 

2 

$ 1.8 
1.6 

1.4 
1.2 

1 

0 0.2 0.4 0.6 0.8 1 





Figure 2. Safety factor profile of the tokamak employed for the calculations of Section 

3. 


In reference [5], Xiao and Catto gave an expression more accurate in the inverse aspect 
ratio expansion. Namely, 

Moo) ^ _ 1 _ 

V7fc(0) 1 + 1.6g%-i/2 + 0.5g2 + 0.36g%V2 ' ^ > 

Both of these results were obtained by using the analytical equilibrium of a large aspect 
ratio circular tokamak which, in our coordinates, is given by 


B = _ — _ ('49') 

l + £cos(27r0)’ ^ ’ 

where Bq is the magnetic field strength at the magnetic axis. The analytical solutions 

(47) and (48) are plotted in figure 1, together with the numerical evaluation of (44) with 

CAS3D-K, for k±pti ^ 1 and a Maxwellian initial condition. We use an axisymmetric 

tokamak with major radius R = 1.7 m, minor radius a = 0.4 m, and q prohle given 

in figure 2. For the CAS3D-K computations, the equilibrium is obtained with VMEC 

employing the aspect ratio and safety factor values just mentioned. The wavenumber 

used in the CAS3D-K calculation is kp = 0.5 and the dimensionless quantity {k±pu)^ 

ranges from 0.0015 in the innermost radial position to 0.0068 in the outermost one. We 

have checked that the residual zonal flow value obtained with CAS3D-K and shown in 

figure 1 does not change if {k±pti)^ is further decreased. The regions of figure 1 where 

the curves agree and where the curves differ are as expected (see the remarks in [32] 

about hgure 3(a) in that reference). The analytical equilibrium of a large aspect ratio 

circular tokamak, used in deriving the equations (47) and (48), differs less from the 

numerical equilibrium obtained with VMEC in radial positions closer to the center. We 

will see in Section 4 that the CAS3D-K results coincide with gyrokinetic simulations of 

zonal flow evolution, in which VMEC equilibria are also used. 
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e 

Figure 3. Magnetic field strength along a field line of the analytical large aspect 
ratio circular tokamak equilibrium, given by equation (49), with e = 0.2, and for the 
numerical equilibrium obtained with VMEC (with R = 1.7 m, a = 0.4 m) at ijj = 0.7. 


As explained above, Xiao and Catto (X-C) also addressed in references [5, 6] the 
extension of the calculation in [3] to short wavelengths. They gave the result 


y’fc(oo) __ ft {1 ~ _ 

Afc(O) |l — Jos^ 


(50) 


for the residual zonal flow in a tokamak at arbitrary wavelengths and with kinetic 
electrons. The expression provided by X-C for the case of adiabatic electrons is 


^kjoo) _ _ ft ~ _ 


(51) 


In order to avoid any confusion, we have to point out that (50) and (51) were not 
derived as the solution of the initial value problem explained in Section 2, but assuming 
that the quasineutrality equation is forced with a source term. The argument of X-C 
can be streamlined as follows. Go back to (41) for the tokamak case (that is, dJJ = 0 
for all particles). X-C consider that hnite orbit width effects do not affect the initial 
condition; i.e. 


S 

J2Zs{Josfs{0)/Fso},= f 

s \ s / ^ 


( 52 ) 
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Figure 4. Comparison of the result in references [5, 6] and the evaluation of (50) with 
CAS3D-K. The parameters of the tokamak are the same as in figure 1. 


Hence, in this approximation, (41) gives 


tpfloo) ^ 


l^s Ts 


( -^ cJJ=0 

{ 1 }. - 


(53) 


The qnasinentrality eqnation at t = 0, (42), can be trivially rewritten as 

Josfs(0)d^v\ . (54) 

/ 

From the qnotient of (53) and (54), one obtains equation (50). Analogously, one can 
obtain (51) from (44). From these manipulations, it is clear that in the X-C calculation 
the charge perturbation at t = 0 can be viewed as a constant source term in the 
quasineutrality equation (see also the remark after equation (57)). 

References [5] and [6] provided analytical evaluations of the right sides of (50) and 
(51) for simplihed tokamak geometries. Since we can directly evaluate the right sides of 
(50) and (51) with CAS3D-K, we will compare the results as an additional check of our 
numerical tool. 

In [5, 6], the analytical equilibrium of a circular cross section, large aspect ratio 
tokamak with safety factor q = 2 and inverse aspect ratio e = 0.2 was used. For the 
calculations with CAS3D-K, we employ the VMEC tokamak equilibrium described above, 
which has similar parameters at '0 = 0.7. At this radial position, the VMEC equilibrium 
satisfies q = 2 and e = 0.2 within an error of 1.5% (e: = 0.197 and q = 2.03). The 


V2fc(0) = 



4 }. 
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Figure 5. Comparison of the result in reference [5] for adiabatic electrons and the 
evaluation of equation (51) with CAS3D-K. The parameters of the tokamak are the 
same as in figure 1. 


difference between the VMEC equilibrium and the analytical one is illustrated in figure 
3, where we compare the magnetic held strength along a held line for both equilibria. 
In the VMEC equilibrium, the value of the magnetic held strength at the magnetic axis 
is Bq = 1.87 T. In general, deviations from circularity are expected in the numerical 
equilibrium because of ehects like the Shafranov shift that are not taken into account in 
the analytical equilibrium. These deviations are smaller for radial positions closer to the 
center. The comparisons for the cases with fully kinetic species (50) and with adiabatic 
electrons (51) are shown in hgures 4 and 5. The agreement is quite good. The fact that 
the curves present some diherences, especially at short wavelengths, is not surprising 
because the equilibria are not identical. As already advanced above, we will see that 
the CAS3D-K calculations agree very well with the results from gyrokinetic simulations 
carried out with Gene coupled to GIST [15] and EUTERPE, that also employ VMEC 
equilibria. It should be noted that further gyrokinetic codes with similar capabilities 
exist. The independently developed code GKV-X [33], for instance, is also able to handle 
VMEC equilibria. 

In the next section, we compare CAS3D-K calculations of the residual zonal flow 
with the results obtained from gyrokinetic simulations. 
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4. Comparison of the residual zonal flow values obtained with CAS3D-K 
and with gyrokinetic simulations 


In this section, we calculate the residual zonal flow as an initial value problem for a 
wide range of radial wavelengths in tokamak and stellarator geometries. We use the 
numerical techniques explained in Section 3 to evaluate the required averages using 
the code CAS3D-K. These calculations will be compared with the results from two 
gyrokinetic codes, showing that the extension of CAS3D-K is faster. 

The residual level for the initial value problem is given by (41), once an initial 
condition /^(O) has been specifled. This initial condition has to satisfy (42). An initial 
condition fulfllling this equation is 

= (55) 

J- Os 

Using (55), we And that the expression for the residual level is 


t^fc(oo) 

(^fc(O) 



(56) 


Similarly, from (44) and (55) we obtain the residual level when using the approximation 
of adiabatic electrons. Namely, 


y’fc(oo) 



(1 - Tos)^ 

Ts 

{ 1 }^- 



(57) 


As explained above, in our numerical evaluations of (56) and (57) with CAS3D-K, we 
assume that in stellarator geometry the trapped trajectories have aU 7 ^ 0 . 

To be precise, the Xiao and Catto formulas (50) and (51) can be obtained from 
an initial value problem calculation by choosing an initial condition /s(0) different from 
ours. However, the initial conditions that recover the X-C results necessarily have 
increasingly fast oscillations along the orbit for increasing k^, and seem of limited interest 
for the analysis of turbulence simulations. For this reason, we choose a different initial 
condition that is in our opinion more relevant. 

The comparison with gyrokinetic simulations will be carried out by using the 
Gene and EUTERPE codes. Gene [12, 13, 14, 15] is a Eulerian gyrokinetic 5f 
code which can be run in radially global, full flux surface or flux tube simulation 
domains. The code can use adiabatic or kinetic electrons and is able to deal with 
tokamak and stellarator geometries. In the Gene simulations, we calculate the zonal 
flow response for a wide range of radial wavelengths, using both adiabatic and kinetic 
electrons. EUTERPE [16, 17] is a global 6f gyrokinetic code in 3D geometry with 
a Lagrangian Particle In Cell (PIC) scheme. In the simulations with EUTERPE, a 
k±Pts < 1 approximation is employed in the quasineutrality equation that limits the 
range of wavelengths for which we can carry out the calculations. With EUTERPE, 
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Figure 6. Safety factor profiles employed in Section 4 for the tokamak (solid line) and 
W7-X (dashed line) calculations. 


we have been able to simulate adiabatic electrons and also kinetic heavy electrons. 
All the gyrokinetic simulations shown in this work are linear and collisionless, with a 
plasma that, unless stated otherwise, consists of only two species: singly charged ions 
and electrons, s = {i,e}. We take flat density and temperature profiles with the same 
values for both species. 

In EUTERPE an initial condition proportional to sin(/c^'0) is used for the perturbed 
distribution function. After the first time step, a zonal perturbation to the potential with 
the same radial dependence (p(0) oc sin(/c^'^) appears, that is used as the initial zonal 
flow. For the implementation of the initial condition in EUTERPE, the Bessel functions 
Jqs and Fos are approximated to lowest order in k±pts <C 1. The initial condition in 
EUTERPE is then 

Fsi{0) = sinikpi;) F,o. (58) 

Since Gene works in Fourier space for the radial coordinate, we initialize the perturbed 
distribution function with only one radial mode which produces a potential with a single 
mode of unit amplitude. 

4.I. Tokamak 

First, we compare gyrokinetic simulations and CAS3D-K calculations in tokamak 
geometry. We use an axisymmetric device with major radius R = 0.95 m, minor 
radius a = 0.25 m, and q profile given in figure 6, whose equilibrium is determined 
by VMEC. We use flat temperature profiles with Tj = Tg. The residual levels obtained 
with EUTERPE, CAS3D-K and the flux tube version of Gene are shown in figure 7 for 
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Figure 7. Residual zonal flow for the initial value problem in an axisymmetric large 
aspect ratio tokamak with major radius R = 0.95 m, minor radius a = 0.25 m and q 
profile given in figure 6. The values predicted by R-H and X-C (equations (47) and 
(48), respectively) are also shown for comparison. 


a radial position -0 = 0.25. We show the calculations with fully kinetic species and also 
using the approximation of adiabatic electrons. The results of the gyrokinetic codes 
have been obtained htting the temporal evolution of the potential to an exponential 
decay model 

</9(t)/(y9(0) = i? + 74exp(—^t). (59) 

The results with CAS3D-K correspond to the evaluation of the equations (56) and (57). 
From hgure 7, we can see that the agreement among the results of CAS3D-K, EUTERPE 
and Gene is excellent. When evaluating the residual level with gyrokinetic codes, a 
certain variability in the results must be assumed. This variability comes from the 
htting method (smaller than 1%), the discretization in phase space and the control of 
the numerical noise, among other factors. All the results in this work obtained with 
Gene show variations smaller than 10%. In any case, hgure 7 shows that the residual 
level obtained by the three independent methods coincides within a margin smaller than 
this quantity, which gives us conhdence to consider that the overall error is quite small. 

As can be seen in hgure 7, the residual value has local maxima centered at the 
scales of the electron and ion Larmor radii. In the long-wavelength limit, k^pu -C 1, the 
residual level in a tokamak does not depend on /cy. Its value is well predicted by (47), 
and even more accurately by (48), for a large aspect ratio tokamak with circular cross 
section. For the VMEC equilibrium used here, these predictions (also indicated in hgure 
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Figure 8. The same evaluation with CAS3D-K of equation (57) as in figure 7, with 
the adiabatic electron approximation, and the evaluation of equation (56) with fully 
kinetic species for different values of r. 


7) are not so accurate for the reasons discussed in the paragraph below equation (48). 
At very short wavelengths, k_\_pte > 2, the residual value approaches zero as 
when using fully kinetic species and also for the adiabatic electron approximation. In 
figure 7, it is shown that the adiabatic electron approximation in tokamaks is good for 
k_LPte ^ 0.1. In hgure 8, we reproduce the results with CAS3D-K in hgure 7 together 
with the evaluation of (56) for different values of r := Tg/Tj. On the electron scale, 
the results with the adiabatic electron approximation and with kinetic electrons only 
coincide in the limit r !:§> 1. Due to the reasons pointed out above, the simulations 
with EUTERPE have been carried out only for k±pti < 1. Finally, it is obvious that the 
results of the forced system of figures 4 and 5 and the initial value problem of hgure 7 
behave in a completely different way for k_i_pti > 1. 

4-2. Stellarator 

Now, we turn to stellarator geometry. We use an equilibrium for the standard 
conhguration of the stellarator W7-X obtained with VMEC. The q prohle is given 
in hgure 6 and we take hat density and temperature prohles with T* = Tg. In 
hgure 9, calculations of the residual level with CAS3D-K, EUTERPE and the full hux 
surface version of Gene are shown for 'ip = 0.25. Two curves correspond to CAS3D- 
K computations, one using adiabatic electrons (see equation (57)) and the other one 
using kinetic electrons (see equation (56)). In hgure 9, the results of the gyrokinetic 
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Figure 9. Residual zonal flow level for the initial value problem in the standard 
configuration of the W7-X stellarator at '0 = 0.25. The q profile is shown in figure 6. 


simulations were obtained employing both the approximation of adiabatic electrons and 
fully kinetic species with Gene, whereas only calculations with adiabatic electrons are 
shown for EUTERPE. These results have been htted to an exponential decay model (59) 
to get the residual value. Similar results can be obtained with an algebraic decay model 
as suggested in reference [10]. The results of CAS3D-K show remarkable agreement with 
both gyrokinetic codes. 

As explained and quantihed at the end of this section, the gyrokinetic simulations 
with kinetic species are much more demanding in terms of computational resources than 
those with adiabatic electrons. Global simulations with EUTERPE using fully kinetic 
electrons in stellarator geometry would require an extremely large computing time. This 
time can be reduced by increasing the mass of the species involved. We have calculated 
for deuterium ions and kinetic heavy electrons s = with tue = 400 me and 

Td = Te- The results are shown in figure 10 where we compare the residual level 
calculated with EUTERPE and CAS3D-K aX, ip = 0.25. The results with adiabatic 
electrons shown in this figure are exactly the same as those with adiabatic electrons in 
hgure 9, obtained for hydrogen ions. Note that, with adiabatic electrons, as the residual 
level only depends on {k±pti)^, the curves for hydrogen or deuterium ions are exactly 
the same. 

In figure 9, like in tokamaks, we find local maxima of the residual level centered 
around the scales of the electron and ion Larmor radii. However, at long wavelengths, 
k±pti 1, the residual level as a function of k_i_pti behaves very differently in tokamaks 
and in stellarators (see, for example, hgures 7 and 9). This can be easily understood 
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Figure 10. Residual zonal flow level for the initial value problem in the standard 
configuration of the W7-X stellarator a,t p = 0.25, with deuterium ions {D) and kinetic 
heavy electrons [E) (jnE = lOOrrie) and also using the approximation of adiabatic 
electrons. 
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Figure 11. The same evaluation with CAS3D-K of equation (57) as in figure 9, with 
the adiabatic electron approximation, and the evaluation of equation (56) with fully 
kinetic species for different values of r. 












Residual zonal flows in tokamaks and stellarators at arbitrary wavelengths 


22 


9- 


9- 


0.01 0.1 1 



Figure 12. Range of validity of the long-wavelength (LW) approximations, when 
using fully kinetic species (60) and with the approximation of adiabatic electrons (61), 
compared to the exact expressions (56) and (57), in the standard configuration of the 
W7-X stellarator at "0 = 0.25 and with = Tg. 


by expanding (56) and (57) in kj^pu -C 1. The nnmerator of these expressions scales 
qnadratically with k±pti in tokamaks and stellarators. The difference comes from the 
denominator. In a stellarator, the denominator is non-zero when k±pti = 0. However, in 
a tokamak the denominator scales qnadratically with k^pu. The denominator has been 
often related to the shielding effects of collisionless classical and neoclassical polarization 
cnrrents [3, 4, 5, 6]. 

It is worth giving explicitly the k^pu 1 expansions of (56) and (57) in a stellarator 
and discnssing a stellarator specific point in detail. The lowest order term of (56) gives 


ipk{oc) (1 - et) 4 

^ (1 + r./(z|r.)) + 


MO) f,- . 

where et = is the fraction of trapped particles. Here, the snperindex 

“trapped” means that the phase space integration is performed only over the trapped 
region. However, if we nse the approximation of adiabatic electrons, (57), we hnd 


iffloo) _ (1 ~ ^i) 


+ 0{{k±pti)^^). 


(61) 


Vkifl) et 

Hence, in stellarators, the adiabatic electron approximation gives and incorrect residnal 
zonal flow, even at long wavelengths. This has been pointed ont in references [7, 8] 
and is confirmed by the calcnlations shown in fignres 9 and 10. The cnrves in figure 
11 for different values of r quantify the error of the adiabatic electron approximation 
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for any wavelength. As can be seen in this hgnre, the residnal level obtained with this 
approximation, only coincides with that obtained with fnlly kinetic species in the limit 
r 3> 1. In hgnre 12, we plot the cnrves in hgnre 9 corresponding to CAS3D-K together 
with the evalnation of their expansions to lowest order in k_i_pti 1, (60) and (61). 
It is clear from hgnre 12 that (60) and (61) are good approximations of (56) and (57) 
respectively for k_i_pu < 0.2. 

We point ont that at scales comparable to the ion Larmor radins, k^pu ~ 1, the 
residnal level appears to be larger in stellarators than in tokamaks (see, for example, 
hgures 8 and 11). In order to discard trivial explanations, we have stndied with CAS3D- 
K the residnal level in a tokamak conhgnration with the same aspect ratio and q prohle 
as those of the standard conhgnration of W7-X and the resnlts are mnch closer to 
the tokamak case than to the stellarator resnlts. This is not snrprising becanse, as 
shown in reference [6], not only the aspect ratio bnt also shaping ehects like elongation, 
triangnlarity and Shafranov shift, among others, ahect the residnal level. We leave for 
a fntnre work a detailed stndy of the magnetic conhgnration inhnence on the residnal 
level. 

4 . 3 . Simulation conditions and computational time 

The converged resnlts shown in this work reqnire disparate compnting resonrces, 
depending on the code and the physical problem. The relevant nnmerical parameters 
and the compntational resonrces reqnired by each code are described below. 

In CAS3D-K, we nsed 256 points for the integration over the velocity coordinate v, 
with 0 < n < 4:7rvts- For the integration over the A coordinate, we nsed 72 integration 
points in the tokamak and 24 in W7-X. Along the held line, we nsed 32 points at long 
wavelengths and np to 4096 for short wavelengths. This resolntion allows the correct 
integration of the highly oscillatory fnnctions. Thanks to axisymmetry, in a tokamak 
all held lines on a hnx surface are equivalent. In W7-X, we used 1024 held lines to cover 
the hnx surface for passing particles. For the evaluation of 6s in the stellarator case, all 
the modes with |m| < 8 and |n| < 8 were retained. The calculations were carried out 
in the EULER cluster at CIEMAT, equipped with Xeon 5450 quad core processors at 3 
GHz and 4XDDR Inhniband network. 

In Gene, a ID spatial grid along the held line {z coordinate) is used in the tokamak 
cases while in stellarator simulations a 2D spatial grid in coordinates {y, z) is used to 
describe a full hnx surface {y is the coordinate along the binormal direction). In velocity 
space, a 2D grid in parallel velocity and magnetic moment coordinates {v\\,p) is used in 
both the tokamak and the stellarator cases. The resolution of the spatial and velocity 
grids used are given in table 1 together with the time step and the total simulation 
time for each case. Times are given in l/kla units, with VLg = a/vte^ where we recall 
that a is the minor radius and Vte is the thermal velocity of electrons. The GENE 
simulations were run in HYDRA [34], equipped with Intel Ivybridge at 2.8 GHz and 
SandyBridge-EP at 2.6 GHz processors interconnected by Inhniband FDR14. 
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Table 1. Numerical parameters used in GENE simulations. The time step {^tc) and 
the total simulation time {Tq) are given in units, with tic = ajvte- 



long wavelength {k±pu < 1) 

short wavelength {k±pu > 1) 


adiabatic e 

kinetic e 

adiabatic e 

kinetic e 

nz 

64 

64 

128 

256 


128 

1024 

1024 - 2048 

2048 - 4096 


40 

40 

40 

40 tokamak 

Ate 

1-6 

0.03-0.06 

0.2 - 1 

0.02 

Tg 

4000 - 15000 

3000 - 10000 

4000 - 10000 

150 

Uy 

64 

64 

64 

64 

Hz 

256 

128 

128 

128 

^^11 

128 

256 - 512 

128 

128 — 512 stellarator 

rif. 

20 

20 

20 

10-20 

Ate 

4-8 

0.06 

0.2-5 

0.06 

Tg 

100000 - 200000 

40000 

300 - 5000 

10000 - 550000 


Table 2. Numerical parameters used in EUTERPE simulations. The time step {AIe) 
and the total simulation time (Te) are given in units, with Qe = eB* jm. ^This 
range corresponds to calculations with deuterium ions and heavy electrons. ^^Only for 
the shortest-wavelength case. 


long wavelength 

tokamak 

stellarator 

{k±pu < 1) 

adiabatic e“ kinetic e“ 

adiabatic e 

kinetic e f 

nsE 

32 - 192 — 

64 - 192 

32-96 

^Oe ^ ^4>e 

16 X 16 — 

32 X 32 

16 X 16 

of markers 

40 M - 240 M — 

40 M - 240 M 

40 M - 120 M 

hlter cutoh 

5 — 

5,10b 

5 

AtE 

10-5 — 

50-20 

0.5 

Te 

60000 — 

400000 

45000 


In EUTERPE, the electric potential is represented in a 3D spatial grid in PEST 
coordinates {se,0e,(I^e) whose radial resolution must be large enough to correctly 
represent the potential perturbation. The number of markers was set according to the 
grid resolution to maintain the ratio of markers per grid cell approximately constant. A 
low-pass squared hlter in Fourier space k^pj^) is used to reduce the noise. In table 
2 the resolution of the spatial grid ), the number of markers, the hlter 

cutoff, the time step and the total simulation time used for each case are given. Times 
are given in I/VLe units, where VLe = eB*/rui, e is the elementary charge, ruj is the ion 
mass and B* is the average of the magnetic held along the magnetic axis. The EUTERPE 
simulations were carried out in EULER and MareNostrum III [35], equipped with Intel 
SandyBridge-EP processors at 2.6 GHz and Inhniband FDRIO interconnection. 

In table 3 we illustrate the computational cost, in total CPU core hours (that is. 
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Table 3. Estimated CPU time (total core hours) to obtain the residual zonal flow 
value with the different codes. We give estimations for tokamaks and stellarators; 
for adiabatic and kinetic electrons; and for long and short wavelengths. ^This range 
corresponds to calculations with deuterium ions and heavy electrons. UPor very small 
wavenumbers, Gene computes the residual zonal flow in a tokamak in approximately 
0.5 CPU hours. 


long wavelength {k^pu < 1) short wavelength {k^pu > 1) 



adiabatic e 

kinetic e 

adiabatic e 

kinetic e 


CAS3D-K 

1.5 

3 

15 

30 


Gene 

lOtt 

1300 

150 

200 

tokamak 

EUTERPE 

7000 

— 

— 

— 


CAS3D-K 

0.5 

1 

4 

8 


Gene 

2000 

40000 

200 

250000 

stellarator 

EUTERPE 

20000 

800001 

— 

— 



the time summed up over all the cores employed in the simulation), for the different 
codes and cases studied. Of course, the values shown in table 3 are simply indicative, as 
they depend on the numerical details of the simulations and the type of CPU employed 
in each calculation. In addition, a systematic analysis of the optimal resolution to carry 
out the computations with each code has not been performed. The main conclusion that 
we can extract from table 3 is that determining the residual zonal flow with CAS3D-K 
is faster than with Gene and EUTERPE. This is specially true for stellarators. The 
reason is that in stellarators only passing particles contribute to the residual value, while 
in tokamaks also trapped particles count, and trapped trajectories typically demand a 
more careful numerical treatment than passing ones. Whereas the CAS3D-K calculation 
simply drops the contribution from the trapped region, the gyrokinetic runs simulate 
all trajectories. However, we can see from table 3 that the computational cost when 
using the gyrokinetic codes is higher in stellarator geometry because it requires increased 
resolution in phase space to obtain converged results. 

We observe that EUTERPE, a 3D global code, requires much more CPU time than 
Gene, particularly in the tokamak case, in which the flux tube version of Gene is 
used. The reason is that EUTERPE simulates the whole plasma while Gene is here 
operated in a radially local limit. The computational cost with EUTERPE increases 
with k±pti, because more flux surfaces have to be considered as k± increases to properly 
resolve the radial structure of the potential in all the plasma volume. In EUTERPE the 
different values of kj_pu at a given radial position are obtained by keeping the value of 
Pti (determined by the ion mass, the temperature and the magnetic field) and varying 
the value of k_t^. In CAS3D-K, the resolution at short wavelengths must be increased to 
correctly calculate the highly oscillatory functions related to the finite orbit width and 
the finite Larmor radius effects. 

The analytical expression obtained by Rosenbluth and Hinton in reference [3] , given 
by equation (47), has been largely used as a linear benchmark for gyrokinetic codes in 
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tokamak geometry and in the long-wavelength limit. The results presented in this work 
show that CAS3D-K can be used to perform those benchmarks not only in tokamak 
geometry but also in stellarator geometry and for arbitrary wavelengths. Examples of 
such benchmarks are given in figure 7 for the global code EUTERPE and the flux tube 
version of Gene in tokamak geometry, and in figure 9 for EUTERPE and the full flux 
surface version of Gene in stellarator geometry. 

Finally, it is a matter of fact that the tokamak calculation including a source term 
in the quasineutrality equation [5, 6] has become quite popular in the literature. Just 
for completeness, we give an analogous calculation for the stellarator in Appendix B 
using gyrokinetic simulations and CAS3D-K. We also show the results for the tokamak 
using gyrokinetic codes (these were not included in Section 3). 

5. Conclusions 

In this work we have treated analytically the linear collisionless zonal flow evolution as an 
initial value problem, and derived expressions (see (41) and (44)) for the residual value 
that are valid for arbitrary wavelengths and for tokamak and stellarator geometries. The 
expressions (41) and (44) involve certain averages in phase space, and also the solution 
of magnetic differential equations, that cannot be evaluated analytically except for very 
special situations. We have extended the code CAS3D-K to evaluate such expressions 
in general. We have tested the extension of the code by comparing its results with 
analytical formulae available in the literature for simplified tokamak geometry [5, 6]. 
These tests are given in figures 1, 4 and 5. 

Then, we have computed the residual zonal flow level in tokamak and stellarator 
geometries for a wide range of radial wavelengths, using both the approximation of 
adiabatic electrons and fully kinetic electrons. We have compared the results of CAS3D- 
K with those obtained from two gyrokinetic codes: the global code EUTERPE and the 
radially local versions of Gene (full flux surface and flux tube). The comparisons are 
shown in figures 7, 9 and 10. 

A stellarator specific effect has been discussed in detail. Namely, the fact that the 
adiabatic electron approximation gives incorrect zonal flow residuals even for k±pti -C 1, 
unlike in tokamaks. This effect has also been confirmed by means of gyrokinetic 
simulations. This is shown in figures 9, 10 and 11. 

Finally, we stress the efficiency of our method to determine the residual zonal flow. 
Gyrokinetic simulations with Gene and EUTERPE to obtain the residual level are 
computationally expensive, especially with fully kinetic species, in the short-wavelength 
region, and in stellarator geometry. On the contrary, the calculations with CAS3D-K 
are less demanding (see table 3), particularly in stellarator geometry. These results 
show that CAS3D-K is a useful tool to calculate fast and accurately the residual level in 
any toroidal geometry and for arbitrary wavelengths. This code is even more useful in 
stellarator geometry as kinetic electrons must be considered to correctly calculate the 
residual level. It can also provide a good benchmark for gyrokinetic codes. 
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Appendix A. Magnetic differential equation 

In this section, we give the solution 5s to the magnetic differential equation (23). We use 
Boozer coordinates. That is, we assume that are such that the contravariant 

form of B is given by (5), and its covariant form is given by 

B = hve - 4VC + fl (V', 0, C) (A.i) 

The square root of the metric determinant is — /pT()/R^. Here, = Rifl’) 

and Ip = Ip{fl>) are the toroidal and poloidal currents, respectively. In Boozer 
coordinates, the radial magnetic drift reads 

[Ip^e + pip. (A.2) 

The parallel gyroradius is dehned as p|p = n||/Bs, whereas Th = B^/{v\fl!'p). 


Appendix A.I. Magnetic differential eguation for passing particles 


For passing particles, Us = 0. The magnetic differential equation to be solved is then 
given by 

n||b ■ = (Us, (A.3) 

which in Boozer coordinates becomes 


{flt'pde + ds — — {Ipde + hdfl P|p- 
This equation is easily solved in Fourier space, giving 


m,n7^0 


mip -|- nit 
m\k' -I- nTl 


(P|b)r 


.^2'K\{m6-\-nQ) 


where the coefficients {p\\s)mn are dehned by 


P\\s — 




^27ri(m0+n(^) 


m^n 


(A.4) 

(A.5) 


(A.6) 
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and C is a constant. By choosing the integration constant as C = —Jp(p||s)oo/'hp, we 
have 


5. = -^ 

% 


J V™ + m 




(A.7) 


This solntion, already given in [9], is valid for any toroidal geometry. For axisymmetric 
tokamaks, it rednces to 5s = —Ip p\\s/^'p- 

Appendix A.2. Magnetic differential equation for trapped particles 

For trapped particles, the magnetic differential eqnation is given by 

n||b ■ VSs = uJs —uTs- (A-8) 

In coordinates {V’, 6*, a}, with a = ( — q9 (and {fj, 9, ((■} being Boozer coordinates), one 
has 


ta., = 




f [Ipde + {It (lff)9a] P\\si 


and eqnation (A.8) then reads 

'^b ^9 A ~ 


where 






p 


n ^daP\\s - n ^dap\\s 


(A.9) 


(A.IO) 


(A.ll) 


The coordinate along the magnetic field line, 9, is not monotonic over the periodic orbit 
delimited by the bonnce points 9bi and 9b2- We define a monotonic coordinate r by 


when cr > 0 


T : = 


C \^b\d9' 

Tb/2 — \Tb\ d9' when a < 0, 


with 


Tb = 2 \Tb\ d9. 

J 6»i,i 

Then, the solntion of (A.IO) can be easily written as 

= -^Ph + I 2 „dT'. 


(A. 12) 


(A. 13) 


(A. 14) 


In order to give a more explicit expression for (A. 14), we use periodicity in r and write 

(l^sa); (A.15) 

with Ub := 271/% and 

/■h,/2 

/ C(;<jQ,(r) cos (/a5;,r) dr. (A.16) 

'o 


{^sa)i — {Tb/2d) 
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Here, the fact that Usai^) is even in r has been employed. Finally, 



(A.17) 


Note that for axisymmetric tokamaks (A.17) simply gives 5s = —IpP\\sl^'p- 

Appendix B. Residual zonal flow with a source term in the quasineutrality 
equation 

Figures B1 and B2 are analogous to figures 4 and 5, but this time we employ the 
equilibrium and parameters of Section 4, and show the results obtained by both CAS3D- 
K and Gene. The gyrokinetic simulations have been carried out by taking vanishing 
initial condition and adding a constant source term to the quasineutrality equation. 

A formulation of the residual zonal flow problem similar to that given in [5, 6] gives, 
in the stellarator case. 


7’fc(oo) 

V5fc(0) 



(B.l) 
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Figure Bl. The evaluation of equation (50) with CAS3D-K for the tokamak of Section 
A at p = 0.25 is shown. The corresponding simulation with Gene including a source 
term in the quasineutrality equation is also plotted. 
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Figure B2. The same calculations as in figure Bl, but employing adiabatic electrons. 
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Figure B3. The results of the forced case for the standard configuration of the 
stellarator W7-X at tp = 0.25. 


if all species are kinetic. In the approximation of adiabatic electrons, one has 


g^k[oo) 


Ts 


{1}^ - 


Os 


1 ^ 5=0 


(B.2) 
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The evaluation of these expressions with CAS3D-K, for the standard conhguration of 
the stellarator W7-X and the parameters detailed in Section 4, is shown in B3. The 
results for Gene simulations are also plotted. 
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